iDesign = 2;

%% test baseline for each subject
for d = 9:length(ALLEEG)
    
    TMPEEG = eeg_checkset(ALLEEG(d), 'loaddata');
    TMPEEG = pop_rmbase(TMPEEG, [-1000 -10]);
    [tmp tbeg] = min(abs(TMPEEG.times+1000));
    [tmp tend] = min(abs(TMPEEG.times+10));
    
    nstim = [0 0 0 0];
    for t = 1:TMPEEG.trials
        datatrial = TMPEEG.data(:,tbeg:tend,t);
        datatrial = mean(datatrial,2);
        
        if 
        switch(nstim)
        
        if 
        
        
        
    % select data
    % -----------
    nChan = 32;
    for iStim = 1:3
       nTrials = size(datStim{iStim}.chan1,2);
       allSubjdata{iStim} = zeros(nChan, nTrials);
       for iChan = 1:nChan
           tmpdata     = getfield(datStim{iStim}, [ 'chan' num2str(iChan) ]);
           tmpdata     = mean(tmpdata(timeIndices,:));
           allSubjdata{iStim}(iChan,:) = tmpdata;
       end;
    end;
    
    % concatenate data for the different stimulus types
    % -------------------------------------------------
    allStimTypes = 
    [ allSubjdata{:} ];    
    
    3:length(STUDY.design(2).setinfo)
    
    EEG1 = ALLEEG( d );
    EEG2 = ALLEEG( d+1 );
    EEG1.data = eeg_getdatact(EEG1, 'verbose', 'off');
    EEG2.data = eeg_getdatact(EEG2, 'verbose', 'off');
    
    naccu = 3000;
    R     = zeros(EEG1.nbchan, naccu);
    GL  = [ [ones(1,EEG1.trials) zeros(1,EEG2.trials)]' [zeros(1,EEG1.trials) ones(1,EEG2.trials)]'  ones(EEG1.trials+EEG2.trials,1) ];
    pGL = pinv(GL);

    range1 = [-2000 -1010];
    EEG1 = pop_rmbase(EEG1, range1);
    EEG2 = pop_rmbase(EEG2, range1);
    range2 = [-1000 -10];

    [tmp t1] = min(abs(EEG1.times-range2(1)));
    [tmp t2] = min(abs(EEG1.times-range2(2)));
    ERPmean  = [ squeeze(mean(EEG1.data(:,t1:t2,:),2)) squeeze(mean(EEG2.data(:,t1:t2,:),2)) ];
    
    % compute regular bootstrap
    
    ERPmean2 = { squeeze(mean(EEG1.data(:,t1:t2,:),2)) squeeze(mean(EEG2.data(:,t1:t2,:),2)) };
    [ t df pvalboot] = statcond( ERPmean2, 'mode', 'bootstrap', 'naccu', 1000);
    
    %ERPmean2(:,1:43)  = ERPmean1(:,1:43)+2;
    %ERPmean2(:,44:94) = ERPmean1(:,44:94)-2;

    % compute R at all time and channels
    % ----------------------------------
    for c = 1:EEG1.nbchan
        for ind = 1:naccu+1
            if ind == naccu+1
                tmp = ERPmean(c,:);
            else
                tmp = shuffle(ERPmean(c,:));
            end;
            params    = pGL*tmp';
            m         = corrcoef(GL*params, tmp);
            R(c,ind)  = m(2); %sqrt(mean((GL*params-tmp).^2));
        end;
    end;

    % find max R for all channels
    % ---------------------------
    [tmp maxE] = max(R(:, naccu+1));
    newR = max(R, [], 1);

    % compute bootstrap GLM
    % ---------------------
    orival       = newR(end);
    [surrog idx] = sort(newR);
    [tmp mx]     = max( idx );        

    len = length(surrog);
    pvals = 1-(mx-0.5)/len;
    pvals = min(pvals, 1-pvals);
    pvals = 2*pvals;
    fprintf('%s, %s, p-values: %2.4f (bootstrap %2.4f)\n', ALLEEG(d).subject, STUDY.datasetinfo(d).group, pvals, min(pvalboot));
end;

%% plot best result for richard
figure;plot(ERPmean(maxE,:));
x1 = 1:EEG1.trials;
x2 = EEG1.trials+1:size(ERPmean,2);
m1 = mean(ERPmean(maxE,x1));
m2 = mean(ERPmean(maxE,x2));
n1 = sum(ERPmean(maxE,x1) > (m1+m2)/2);
n2 = sum(ERPmean(maxE,x2) > (m1+m2)/2);
hold on; 
plot(x1, m1*ones(1,length(x1)), 'k');
plot(x2, m2*ones(1,length(x2)), 'k');
